Helical edge states in multiple topological mass domains 
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The two-dimensional topological insulating phase has been experimentally discovered in HgTe 
quantum wells (QWs). The low-energy physics of two-dimensional topological insulators (TIs) is 
described by the Bernevig-Hughes-Zhang (BHZ) model, where the realization of a topological or 
a normal insulating phase depends on the Dirac mass being negative or positive, respectively. We 
solve the BHZ model for a mass domain configuration, analyzing the effects on the edge modes of a 
finite Dirac mass in the normal insulating region (soft- wall boundary condition). We show that at 
a boundary between a TI and a normal insulator (NI), the Dirac point of the edge states appearing 
at the interface strongly depends on the ratio between the Dirac masses in the two regions. We also 
consider the case of multiple boundaries such as NI/TI/NI, TI/NI/TI and NI/TI/NI/TI. 

PACS numbers: 73.43.-f, 72.25.Dc, 73.20.At, 73.21.Fg 
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I. INTRODUCTION 

Topological insulators (TIs) are time-reversal- 
symmetric materials featuring a topological phase 
characterized by a Z2 topological invariant*!- 2 -. In two- 
dimensions (2D), they exhibit the quantum spin Hall 
(QSH) phased. The QSH phase has been theoretically 
predicted 4 - and experimentally realized in HgTe/CdTe 
QWs^. The crucial ingredient of this narrow gap 
semiconductor material is the inverted band structure of 
HgTe. Similarly, 3D TIs supporting chiral fermions as 
surface states have been proposed and observed^—. 

In HgTe/CdTe QWs, the topological phase is deter- 
mined by the sign of the Dirac mass M. The gap be- 
tween the El (s-like) and the HI (p-like) subbands at 
the r point is given by 2|M|. The only experimentally 
accessible parameter tuning the Dirac mass from normal 
(M > 0) to inverted (M < 0) is the thickness of the HgTe 
QW. In particular, a topological transition from the nor- 
mal to the topological insulating phase takes place when 
the QW thickness is increased above the critical thickness 
tc = 6.3 mr£. Recently, electrically driven topological 
insulating phase transitions have been proposed in het- 
erostructures with gate tunable conduction-valence band 
energy separation. In particular, in Ref. [TU a type-II 
InAs/GaSb/AlSb QW was proposed and recent experi- 
ment o 13 i 14 provided the first evidence pointing towards 
the presence of a topological insulating phase in these 
structures. In Ref. 15, double QW structures composed 
of narrow gap semiconductors are considered featuring a 
tunable topological transition with the application of a 
gate bias of the order of the gap of the individual trivial 
QWs. Both proposals pave the way to 2D systems where 
mass domains are designed with lithographic gates to cre- 
ate topological and normal regions. For example, ring- 
shaped TI regions are particularly interesting for the pe- 
culiar properties of their confined edge states, which can 
be controlled with a threading mag nctic Bxn^ML. On the 
other hand, the Dirac mass term in single HgTe QWs is 



related to the QW thickness and therefore thickness fluc- 
tuations of less than 1 nm height can accidentally deter- 
mine the formation of mass domains alternating TI re- 
gions, where M < 0, to NI regions, where M > 0J£- These 
phenomena can be especially relevant for near zero-gap 
HgTe QWsM. 
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FIG. 1: (Color online) Schemes of multiple topological mass 
domains analyzed in the present paper. According to the 
classification of the topological phase of each region into a TI 
or a NI we have M < or M' > and M" > 0, respectively. 

The change in the topological invariant Z2 between 
two systems determines the presence of ID helical edge 
states running along the boundary between the TI and 
the NI regions, a phenomenon referred to as the bulk- 
boundary correspondence. Such edge states are topolog- 
ically protected against single particle elastic backscat- 
tering (as long as time reversal symmetry [TRS] is pre- 
served) and are particularly interesting for their spin and 
charge transport properties. In the literature, these edge 
states are generally obtained by solving the Hamiltonian 
with hard- wall boundary conditions (BCs), i.e. by im- 
posing that the wave function vanishes at the interface 
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between a TI and a normal medium (an exception is 
the recent proposal of natural boundary conditions^ 1 -). 
While this condition is appropriate to treat the inter- 
face between TIs and the M — > +00 vacuum, this is 
by no means a good approximation when dealing with 
electrically-induced mass domain o 12 ' 15 , where the normal 
regions have a finite positive Dirac mass. Finite mass- 
domains in Dirac systems have first been solved in the 
context of zero energy bound states in the 1+1 Dirac 
equation by Jackiw and Rebbi 2 ^ and for interface states 
in band-inverting contacts based on HgCdTe and PbSnTe 
23 i 24 . More recently, finite mass-domains were proposed 
to induce valley-polarized metallic states in biased bilayer 
graphene 2 ^. 

In the present paper we solve the edge states of a 
HgTe/CdTe QW-based TI for the case of soft-wall BCs, 
appropriate for describing systems with Dirac mass do- 
mains, where the wave function does not vanish at the 
interface, but its continuity and the continuity of its nor- 
mal derivative are instead required. In Section II, we 
briefly review the BHZ model and describe the method 
used to solve the edge states in the topological mass do- 
mains. In Section III, we deal with a single NI/TI inter- 
face schematized in Fig. HJa), with both hard- wall and 
soft-wall BCs. We show that soft-wall BCs quantitatively 
change the dispersion curves [Fig. 2(a)] with respect to 
the case of hard-wall BCs. In particular, the hard-wall 
limit is only reproduced in the limit M' — > 00. Even 
for M' = 10 4 meV, which is of the order of magnitude 
of the electron extraction work-function of a crystal, an 
appreciable deviation from the hard- wall limit is still ob- 
served. The density profile of the edge state bound to 
the mass-domains (Fig. 3) is qualitatively affected by the 
soft-wall BCs, which allow it to extend into both the 
TI and the NI regions. In Section IV, we address the 
effect of a finite bulk inversion asymmetry (BIA) term, 
which introduces off-diagonal matrix elements to the oth- 
erwise block diagonal form of the BHZ model, and eval- 
uate its effects on the bulk dispersion curves and on the 
edge state dispersion (Fig. [5j . The effect of BIA on the 
edge states of a single NI/TI is shown to be tiny, so that 
BIA can be safely neglected in the more complex case 
of multiple NI/TI boundaries. In Section V, we analyze 
a system with two TI/NI interfaces. In particular, we 
investigate the case of a TI strip embedded in a normal 
system, sketched in Fig.QJb), and vice versa: the case of 
a TI system where a strip region with M > is present 
[Fig.QTc)] and analyze for the first time the edge coupling 
through the normal region. In both cases, the overlap of 
edge states bound to different boundaries leads to a fully 
gapped edge mode spectrum (Fig. [7] and Fig. [TT]), with 
a minigap exponentially shrinking with the distance be- 
tween the two interfaces. 

Soft-wall BCs can also be used to couple three or 
more edge states, thus, in Section VI, we analyze the 
edge states for a system with three NI/TI boundaries 
[Fig. Did)]. We study how the properties of this system 
vary with the Dirac mass and geometrical parameters. 



We further argue that this system can describe a helical 
edge state at the sample boundary (vacuum/TI interface) 
in the presence of Dirac mass fluctuations in the TI com- 
position giving rise to mass domains with bubbles having 
M > (i.e. normal character) in the bulk of the sam- 
ple whose edge states could interact with the helical edge 
states at the sample boundary. 



II. THE BHZ MODEL 

The spectrum of a HgTe QW near the T-point is effec- 
tively described in its low energy sector— by the 4-band 
model^ 
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k 2 + fc 2 and a is the vector of Pauli 



where k = \k\ - 

matrices associated with the band-pseudospin degree of 
freedom (band E\ or -Hi)2£. H? is represented in the ba- 
sis \H\+), \Ei-), l-Hi-)}, where the E x states 
(J z — ±1/2) are a mixture of the s-like Tq band with 
the F§ light-hole band, while Hi (J z = ±3/2) is basi- 
cally the T% heavy-hole band. For later use in numeri- 
cal simulations, we quote the following choice of param- 
eters: A = 375 meV nm, B = —1120 meV nm 2 and 
D = —730 meV nm 2 . These parameters follow from 
the 8x8 Kane model 2 ^. Without loss of generality we 
also assume C = 0. The Dirac mass M depends on the 
QW thickness and M < corresponds to the inverted 
(QSH) regime whereas M > corresponds to the normal 
regime. In a first approximation, is block diagonal 
in the spin degree of freedom^, for which we define the 
corresponding vector of Pauli matrices f. As we con- 
sider only systems with TRS, we can restrict ourselves 
to the block h(k). Results can be extended to the other 
Kramers block h*(—k) which is related to h(k) by the 
time reversal operation T = iT y aoK, where K is the op- 
erator of complex conjugation. 

The bulk dispersion curves obtained as the eigenvalues 
of Eq. fT} are described by 



E±(k) = e k ± y/(M - Bk 2 ) 2 +A 2 k 2 . 



(2) 



With a standard choice for the TI parameters (like the 
parameters stated above) the bulk dispersion relation dis- 
plays a conduction band minimum (valence band max- 
imum) at fc = with energy E = M (or E — —M). 
However, depending on the values of the parameters in 
Eq. (U), the bulk dispersion curves can also show a "Mexi- 
can hat" behavior. For a detailed analysis of the behavior 
of the 4-band model as a function of its parameters, see 
Appendix [£] and Ref. [H 
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A. Boundary conditions 

We are interested in obtaining the eigenstates of 
Eq. (TTJ) in real space for a semi-infinite geometry, invari- 
ant under translations along the x-axis. For such a sys- 
tem, k y is no longer a conserved quantity and should 
be replaced by the operator —id y . Compatible with 
a fixed energy E and a real k x , the secular equation 
\h(k) — E\ = provides four fc^-modes: 



tf = k& = -k*-F±y/F*-Q* 



with 



F 

Q 2 



A 2 - 2{BM + DE) 
2{B 2 ~ W) ' 
M 2 - E 2 



B 2 -D 2 



(3) 

(4) 
(5) 



With our choice of parameters, k± has one imaginary and 
one real solution for E within the energy range of the bulk 
bands, whereas in the bandgap (— \M\, \M\) both values 
for k y are imaginary. 

For each fcy-mode, one can write the spinors satisfy- 
ing the Schrodinger equation and corresponding to the 
Kramers blocks r = ±1 as 



(6) 



The ratio between the two components is 

A(rk x *^y) 



-M — E — (D — B)(k% + k 2 ) 



The general solution of the Dirac equation with energy 
E and wave vector k x is therefore given by a linear com- 
bination of the four solutions k y = Xk^ (with A, fx = ±) 
obtained from Eq. ([3]), 

*tU*>v) = E Cx]s/>k.M„Ax>v), ( 8 ) 

where we have introduced the index n to refer, in what 
follows, to the n-th mass domain region. 

Here we discuss the general procedure we use to solve 
the BHZ model in a system composed by N mass do- 
mains with parallel boundaries at y = y n with n = 
0,1, ... N. Inside the n-th mass domain with the con- 
dition y 6 (y„_i,?/„), we consider the Dirac mass term 
(and all other parameters) as constant and a general ex- 
pression for the spinor is given by Eq. ©. At y = y n the 
value of the Dirac mass M changes step-like^. 

Hard-wall BCs at y for the domain n are expressed by 



(9) 



meaning that the edge state cannot extend beyond the 
boundary, being subject to a hard-wall confinement. 



Soft-wall BCs between two consecutive domains n—1 and 
n are instead expressed by the continuity of the spinor 
and its normal derivative 



dy^k 



(n-1) 



(x,y n ) 



(10) 



A system of BCs (either soft or hard) can be always 
expressed in a compact form as 

M^(£)c=0, (11) 

where c is a vector containing all the free coefficients 
characterizing the wave function [Eq. in the n = 
1,...N mass domain. M.k !C (E) can be constructed by 
appropriately using either Eq. ([SJ or Eq. (flU| at each 
one of the boundaries of the system. As in standard 
quantum mechanics, the BCs determine the eigenenergies 
_Efe x through the secular equation 



det[M fcx (£;)] =0. 



(12) 



The corresponding eigenspinors we obtain by solving 
Eq. pT|) for the coefficients c. 



III. AN ISOLATED BOUNDARY 

In this section, we consider a single interface between a 
TI (y > 0) with M < and a NI (y < 0) with M' > [see 
Fig. [Ha)] and calculate the resulting helical edge modes. 
(7) In the limit M' — > +oo, we recover the usual hard- wall 
BCs of Eq.®. For a finite positive M' , soft- wall BCs of 
Eqs. (flUl) are employed. 



A. Hard-wall boundary conditions 

We apply vanishing BCs at y = and search for modes 
with energy lying within the bandgap. Only two of the 
four k v solutions of Eq. ^ with positive imaginary part, 
which we define as k± , are normalizable in the region y > 
and contribute to the edge states. When k± are purely 
imaginary, we define k± — k± . When the parameters are 
such that k± are complex (see in Appendix [A] Eq. IA8[) , 
we instead have k± — ±k±. Using Eqs. (|9|) and (fT2|) . we 
obtain the relation 



R 



(13) 



that imposes a strict relation between the energy E and 
momentum k x . Isolating k x terms and squaring twice, 
we arrive after some algebra at the edge mode dispersion 



curve; 
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FIG. 2: (Color online) (a) Energy dispersion relation of helical 
edge states at the interface between a TI region with Dirac 
mass M — — 10 meV and a NI region of mass M' . Solid lines 
and dashed lines correspond to the two Kramers blocks of the 
Hamiltonian Eq. (1). (b) Value of the Dirac point -Efc x =o as 
a function of M' (solid line). The dashed red line shows the 
hard-wall limit M 1 -> +oo, E kx=0 = -MD/B. In the inset, 
edge states' characteristic decay length in Eq. ([3j ] in 

the NI (full curve) and in the TI (dashed curve) plotted as a 
function of M' at k x = are shown. 



First of all, we observe that after the first squaring of 
Eq. (fT5|) . we lose track of r, therefore only one of the ± 
signs in Eq. (fT4")l is actually a solution of Eq. (|13p for a 
given r. More important, because of the second squaring, 
the solutions in Eq. (fT4|) are not always allowed. For 
the usual parameter choice describing TIs, solutions are 
always admissible with M < O^L. For M > 0, the system 
admits no edge modes. 



B. Soft-wall boundary conditions 

As mentioned previously, soft-wall BCs imply that at 
the interface (y = 0) the spinor and its normal deriva- 
tive are both continuous [see Eq. ([T0|)]. For the M' > 



domain, i.e., y < 0, only the k y modes of Eq. ([3]) with 
negative imaginary part are allowed, while for y > pos- 
sible solutions contain the modes with Im(fcj,) > 0. Using 
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FIG. 3: (Color online) Probability density of the edge state 
for k x — as a function of the distance from the NI/TI in- 
terface (y — 0). The different curves represent the hard- wall 
(solid line) and soft-wall (dashed lines) boundary condition 
cases with M' = 10, 100 meV for a fixed TI mass M = —10 
meV. The probability densities for the two blocks r = ±1 are 
degenerate. 



Eqs. (fTU|) and (fTTj) we obtain the 4x4 matrix 
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and R 



/ 

(15) 
calculated 



where the prime stands for both k± an^ ± v T k± 
in the M' domain. By numerically solving dct [IWe^J — 
0, we determine the energy dispersion relation of the he- 
lical edge states. In Fig. 2(a) we show the energy disper- 
sions i?fc x for a NI/TI interface keeping the TI Dirac mass 
M = —10 meV and varying the NI mass M'. Solid and 
dashed lines correspond to helical edge states of the spin- 
blocks r and — r, respectively. The slope of the curves 
(velocity v x ) is not altered by varying M'. However, as 
shown in Fig. 2(b)| the Dirac point rises with increasing 
M' from -Efc x =o = (its limiting value for M' 0) and 
eventually saturates at the value Ek m =o = [ sce 
Eq. (|14p] for M' — s- oo, reproducing the hard- wall case. 
In the inset of Fig. 2(b) we plot the decay length in the 
TI and NI regions as a function of M'. 

Figure [3] displays the corresponding probability densi- 
ties l^fc^r (x,y) \ 2 of edge states at the Dirac point for 
soft-wall BCs with M' = 10, 100 meV (dashed lines) and 
hard- wall BCs (solid line). In both cases, the wave func- 
tions are strongly peaked closely to the interface y = 
and exponentially decaying away from the interface. The 
characteristic decay length of the edge states are given by 
the inverse of the smaller k y mode calculated from Eq. J3|) 
compatible with their eigenenergies and k x values. 



IV. BIA EFFECTS 

In the present section we address how bulk inversion 
asymmetry (BIA) affects the TI Hamiltonian in Eq. ([T]) 
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and the edge mode solutions for hard-wall boundary con- 
ditions using the model of Ref. |33|. We show that such 
BIA introduces a weak non-linearity in the edge disper- 
sion, especially near the bulk band edges. The position of 
the Dirac point is not affected while the velocity of the 
helical particle is only slightly modified. We therefore 
will not consider the effect of BIA in the sections that 
follow. 
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FIG. 4: (Color online) Bulk-energy spectrum at k = includ- 
ing the BIA-Hamiltonian Eq. (|16p as a function of Dirac mass 
parameter M. It shows an anticrossing with gap 2|A|. 

The block diagonal form of the BHZ model is indeed 
an approximation for the low-energy physics of a HgTc 
QW. The presence of BIA introduces a coupling between 
\E1,±) and \Hl,=p) bands. The leading-order BIA per- 
turbation term is expressed as22 




/ 



Hbia = 



-A 



A 



V-A 



A 



(16) 



It preserves TRS and therefore does not affect the topo- 
logical properties of the BHZ model Eq. (HJ 33 , as long as 
the bulk gap is not closed. Note that Eq. ([T6")) introduces 
an anticrossing at k = as a function of M, see Fig. |U 
Such an anticrossing has been also found in Ref. |3J as 
a function of the QW thickness and the inversion cross- 
ing of the Ei- and iJi-bands shifting to finite fc-values. 
Here we show (see Fig. [S]) explicitly, that the helical edge 
states are still present in the model considered and only 
slightly modified, despite the anticrossing at k = 0. 

Let us first rewrite the BHZ Hamiltonian in Eq. <jTJ) , 
including the BIA term in Eq. (1161) , in the following form 



(17) 



H = e k I 4x4 + AfcS + AA„ - M k k z Y, z , 



where we have introduced two sets of unitary and Her- 
mitian matrices 



= T z a x ; 

A- nr. — Tnr.O'i, 



A,, 



T71 (To, 



A- 



-T z a z 

T z (To 



(18) 



with the property that each set separately obeys Pauli 
commutation rules, while elements from the two sets 



commute. Note that the only matrices which are off- 
diagonal in the Kramers block pseudospin (i.e. contain- 
ing t x or Ty) are A x and A y . We now perform the fol- 
lowing unitary transformation in the A-E space which 
warrants Eq. {jXTJ) block-diagonal: 



U = -j= [-i(CyXx ~ CxZ y )A y + E,A Z ] (19) 

with £ = k/\k\. After the transformation the Hamilto- 
nian in Eq. (|17p acquires the following form 



Hj: = e fc /4 X 4+(-A|*|J 4 x4 + AA,)CE-M fc A,E 



h+(k) 
h*_(-k) 



(20) 



where we have introduced the helicity parameter 77 = ±1, 
and have defined 

h v (k) = s k I 2 x2 + (-A\k\ + r?Aj (( x a x - C y a y ) + M k a z . 

(21) 

The helicity of the energy eigenstates is defined in the 
new basis by ?y = (r 2 ), which in the original basis of 
Eq. IPTl) is equivalent to 77 = — (CE X ( X + E„C«) A„). 
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FIG. 5: (Color online) Bulk bands dispersions including the 
effect of BIA terms, where full and dashed lines distinguish op- 
posite helicities. Edge states for a single hard-wall boundary, 
numerically obtained by solving Eq. (|17p are shown in dotted 
lines. The edge states dispersions for the corresponding sys- 
tem without BIA, obtained with Eq. (|14[l . are also displayed 
for comparison (full narrow line). We used the following pa- 
rameters: M = —2 meV, A = 1 meV. 

The bulk dispersion curves obtained from Eq. (|20[) are 
shown in Fig. [5l Similarly to an electronic system in 
the presence of the Rashba spin-orbit interaction, the 
dispersion curves can be classified through the helicity 77. 
For a given wave vector, the effect of the BIA term is to 
lift the degeneracy of the two spin-blocks. 

In order to solve for the edge states of the system with 
BIA, we need to treat Eq. p7|) in real space [we note that 
Eqs. (TT91 and therefore Eq. (|20p are well defined only in 
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momentum space]. We follow the procedure illustrated 
in Section III A[ applying it to the 4x4 Hamiltonian in 
Eq. (fTT|) . In Fig. El we show the edge states obtained 
for a system with M = —2 mcV and A = 1 meV (dot- 
ted lines) and compare them with the edge states of the 
corresponding system with no BIA (full narrow lines), 
obtained analytically with Eq. (|14l) . These BIA terms do 
not change the position of the Dirac point, but slightly 
change the group velocity close to the Dirac point. Away 
from the Dirac point the edge dispersion in Fig. [5] shows 
a weak non-linear distortion, accentuated near the bulk 
band edges. Projecting Hbia onto the unperturbed edge 
states, it is straightforward to show that the effects of 
Hbia are at least of order A 2 , see Ref. l3ol 
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FIG. 6: (Color online) Energy dispersion of a TI strip 
(NI/TI/NI) calculated with hard- wall boundary conditions. 
The solid red lines correspond to L = 100 nm, blue dashed 
curves to L — 200 nm and black thin lines to L = 1000 nm. 



L = 1000 nm (thin full lines), the overlap is negligible, 
the minigap is exponentially suppressed. The dispersion 
curves are linear, jus t two copies of single-interface edge 
modes shown in Fig. 2(a) for M' — > 00. For L = 200 and 
100 nm, the overlap is instead substantial and the edge 
modes anticross at k x = giving rise to a finite minigap. 




FIG. 7: (Color online) Logarithmic plot of the minigap value 
as a function of the width L of the TI strip for soft-wall bound- 
ary conditions (solid lines) with M' = M" = 10, 50, 100 
and 350 meV. The dashed curve corresponds to the hard- wall 
boundary condition. 



2. Symmetric soft-wall boundary conditions 



V. SYSTEMS WITH TWO BOUNDARIES 

A. A strip of TI 

Here, we consider a NI/TI/NI mass domain shown in 
Fig. Hlb), where Dirac masses are M' > 0, M < and 
M" > 0, respectively. For very large M' and M" the 
use of hard-wall BCs is appropriate. For the entire sub- 
section, the Dirac mass for the TI domain is M = — 10 
meV. 



1. Hard-wall boundary conditions 

The case of a TI strip confined by hard- wall BCs has 
been first analyzed by Zhou et alA We briefly com- 
ment in this section some of their results in order to set 
a benchmark for successive extension to soft-wall BCs. 
The TI strip has two pairs of helical edge states (it is 
not topologically protected) exponentially localized at 
the two boundaries which are separated by the width 
of the TI strip L. The decrease of the TI strip width 
leads to a finite overlap of edge modes belonging to dif- 
ferent interfaces originating a minigap (a full gap in the 
edge mode dispersion curves) as shown in Fig. [5] For 



We now consider M' = M" to be finite, adopt soft-wall 
BCs and solve for the edge modes. The edge dispersion 
curves are qualitatively similar to the hard-wall case in 
Fig. El and an anticrossing behavior is found at k x = 0, 
with the opening of a minigap around the Dirac point of 
the corresponding one-boundary edge modes. We note 
that such an anticrossing point scales in energy with M', 
similarly to the behavior of the Dirac point with soft-wall 
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FIG. 8: (Color online) Electronic density of the edge mode of 
a TI strip of width L — 100 nm at k x — 0. Considered NI 
Dirac masses are M' = M" = 10, 100 meV and 00 (hard-wall 
confinement). 
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M' = 1 meV M' = 1 00 meV M' = 1 000 meV 




I Li I L lJ . II I I I l> . I_J I I I I I L 
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k x (nm" 1 ) 

FIG. 9: (Color online) Edge modes of a TI strip with hybrid 
confinement: M" -> +oo and M' = 10, 100, 1000 meV for 
a TI region of width L = 100 nm. The full and dashed lines 
distinguish r = ±1. 
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FIG. 10: (Color online) Hybrid case - one hard-wall (M" — >■ 
+oo) and one soft-wall (M' = 100 meV) with L = 100 nm. 
The panels (a)-(d) correspond to states labeled by a-d in 
Fig. E 



BCs shown in Fig. 2(b) We plot in Fig. [7] the minigap's 



exponential decay as a function of the TI strip width L. 
The dashed line corresponds to the hard-wall case, and 
the soft-wall cases with finite Dirac masses of 350, 100 , 
50 and 10 meV are also shown. With soft-wall BCs, the 
minigap is smaller and subject to have a faster decay, 
whose origin can be easily understood by looking at the 
edge states profile in Fig. [8j Here, we plot the profile 
of the electronic density of the k x — edge modes in 
the direction perpendicular to the boundaries, for L = 
100 nm and for different NI masses ranging from 10 meV 
to oo (hard- wall case). As one increases M' the edge 
modes become more strongly confined into the TI strip, 
thus leading to an enhancement of the overlap between 
the edge modes at the two interfaces (signaled by the 
increase of the probability at y = 0). 



B. Laterally coupled TI edge states 



We consider also the possibility to laterally couple edge 
modes in a TI/NI/TI mass domain, where the overlap of 
the edge modes takes place in the central NI region [see 
Fig- He)]. For concreteness, we assume TI regions with 
equal Dirac masses of M = —10 meV, while the NI mass 
is M' > 0. This situation leads to qualitatively similar 
edge modes as in a TI strip with soft- wall confinement, 
with the opening of a minigap at k x = as shown in 
Fig lllf a). However, from a quantitative point of view, 
edge modes decay differently in the NI region, according 
to their characteristic penetration lengths [see the inset 
in Fig. 2(b)| . In Fig.fTTTb). we display the minigap value 
which exponentially shrinks as a function of the the NI 
mass M'. 



3. Asymmetric confinement 

We consider now NI/TI/NI mass domains which differ 
from the previous ones for the lack of the mirror symme- 
try at the center of the TI strip. This is the case whenever 
A f 7^ M" . In particular we will focus on the prototypical 
situation of a hybrid confinement, when one interface is 
treated with hard- wall BCs and the other with soft- wall 
BCs. In Fig. [9l we present the dispersion curve for a TI 
strip with hybrid BCs, where M" = oo and for M' = 10, 
100 and 1000 meV. Due to the lack of the mirror symme- 
try the edge dispersion curves display minima at finite 
k x . the effect of increasing M' is to partially compensate 
the imbalance between the two boundaries [M" = oo). 
As a consequence the minima shift towards k x — and 
the center of the gap tends to the hard-wall Dirac point 
value in Fig. [2(b)] The gap is also increasing because of 
the stricter confinement from M' . In Fig.[TUl we plot the 
edge modes profile for the case M' — 100 meV of Fig. |H1 
corresponding to the points a, b, c and d. 
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FIG. 11: (Color online) In (a), edge modes for a TI/NI/TI 
mass domain with width of the central NI region of L = 
100 nm and NI mass of M' — 10 and 20 meV are shown. 
In (b), we present the value of the minigap as a function of 
M'. 
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FIG. 12: Edge mode dispersion curves for a NI/TT/NI/TI 
mass domain with three boundaries [see Fig.[TJd)], calculated 
for a width of the first Tl region L = 100, 70 and 50 nm. The 
width of the second NI region is d = 70 nm with Dirac mass 
M' — 10 meV, while the first NI/TI interface is treated with 
hard-wall BCs. The second Tl region is semi-infinite. Full 
and dashed lines distinguish r = ±1. 



VI. THREE MASS DOMAIN SYSTEM 

In this section we analyze the edge states for a mass do- 
main NI/TI/NI/TI with three boundaries [see Fig.QJd)] 
with Dirac masses M" > 0, M, M' > and M, respec- 
tively. For simplicity we keep the same Dirac mass for 
the Tl regions M — —10 meV. The three-boundary sys- 
tem can only be realized if M' is finite, while M" can be 
either finite or oo leading to soft-wall or hard-wall BCs 
at the first NI/TI interface. The edge dispersion curves 
are quite complex and are more easily understood by first 
considering the system assuming uncoupled edge states 
at each of the TI/NI interfaces. Without coupling, at 
each TI/NI interface we expect linear edge modes similar 
to that in Fig. 2(a)| where the energy value of the Dirac 
point depends on the difference of the absolute value of 



the masses between NI and Tl as described by Fig. 2(b) 
When edge states belonging to different TI/NI interfaces 
overlap a minigap is formed due to anticrossing of the 
dispersions. The anticrossing takes place at k x = if the 
two edge modes have equal Dirac point values (i.e. the 
two TI/NI interfaces share the same parameters), other- 
wise the anticrossing happens at finite a k x . 

We focus our analysis on a system where the first NI 
has a very large bandgap (e.g. the vacuum) and send 
M" — > oo (hard- wall BCs). We define L as the width of 
the first Tl domain (the second Tl domain is considered 
semi-infinite) and d the width of the second NI region 
(the one with mass M'). We note that this situation 
is qualitatively analogous to that of a HgTe QW in a 
Tl phase. The first NI/TI interface is the physical edge 
at the interface vacuum/HgTe QW, correctly described 
with hard-wall BCs. The second NI region can be due to 
a large-scale (tens of nanometers in the 2D-plane) fluc- 
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FIG. 13: (Color online) Electronic density of edge states for 
a NI/TI/NI/TI mass domain [see Fig. \V[.d)] with hard-wall 
BCs at the first NI/TI interface (y = 0). The width of the 
first (from left to right) TI domain is L — 70 nm, while that 
of the second NI domain is d = 70 nm. The Dirac mass of the 
second NI domain is M' = 10 meV. The density profiles in 
the panels (a)-(h) correspond to states labeled by letters a-h 
in Fig. [TJ 



tuation in the QW thickness, leading to the appearance 
of a topologically trivial region. HgTe QWs are typically 
grown to have a thickness around the critical value of 
6.3 nm and it is generally sufficient to have a variation 
of the thickness of the order of fractions of a nanometer 
to induce a band inversion into a NI system. 

In Fig. [121 we plot the edge modes for a mass domain 
with M' = 10 meV and <i = 70 nm and three differ- 
ent values of L = 100, 70 and 50 nm. The edge mode 
of the first NI/TI interface would have by itself a Dirac 
point at around 6.5 meV, however the overlap with the 
edge modes from the second TI/NI interface, which orig- 
inates an anticrossing at k x m 0.01 nm -1 , pushes it up 
to around 8 meV. This analysis is confirmed by the fact 
that the electronic density calculated at point a and / of 
Fig. HU for the case with L = 70 nm [shown in Fig. fTSTa) 
and (f)] is strongly peaked near y — 0, corresponding to 
the first NI/TI boundary with hard-wall BCs. The an- 
ticrossing and the minigap opening at k x ~ 0.01 nm -1 
are due to the overlap of edge states between the first 
and the second boundaries as observed in the spinors in 
Fig. nUT b) and (e) . The gap opening at k x = around 
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FIG. 14: (Color online) Edge mode dispersion curves for a 
NI/TI/NI/TI mass domain with three boundaries, calculated 
for a width of the second NI region d = 100, 50 and 20 nm for 
a fixed Dirac mass M' = 10 meV. The width of the first TI 
region is L = 70 nm, while the first NI/TI interface is treated 
with hard-wall BCs, and the second TI region is semi-infinite. 
Full and dashed lines distinguish r = ±1. 
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FIG. 15: (Color online) Edge mode dispersion curves for a 
NI/TI/NI/TI mass domain with three boundaries calculated 
for a Dirac mass of the second NI region M' — 10, 5 and 
20 meV. The width of the first TI region is L = 70 nm, the 
width of the second NI region is d = 70 nm, while the first 
NI/TI interface is treated with hard- wall BCs, and the second 
TI region is semi-infinite. Full and dashed lines distinguish 
r = ±1. 



E = is due to the overlap of the edge modes of the 
second and third boundary and the edge states resemble 
the case of the lateral coupling of edge states through a 
narrow NI region with mass M' — 10 meV (analyzed in 
Fig. [Til and related text), as can be seen in Fig.[T3ld) and 
(g) . Spinors in Fig. [T5tc) and (h) , which belong to the 
points c and h of Fig. [T^] for the case with L = 70 nm, 
resemble the edge modes of a symmetric TI strip with 
soft-wall BCs away from the minigap region of the dis- 
persion curve. 

The effect of the overlap of edge modes belonging to 
the second and third boundaries can be analyzed by vary- 
ing d, which is done in Fig.[]3] Decreasing d from 100 nm 
to 50 nm accentuates the anticrossing behavior at k x = 
around E — meV due to the overlap of edge states 
belonging to the second and the third interfaces. Other 
features of the edge dispersion are only slightly affected. 
When the second NI region is thinner (see case with 
d = 20 nm) the overlap of the edge states bound to it is 
so strong that they are energetically pushed into the bulk 
spectral range. As a result, the helical edge modes be- 
longing to the first NI/TI interface (Dirac point at around 
6.5 meV) are hardly affected by the presence of a second 
thin NI region. 

The effect of varying the NI mass M' is instead shown 
in Fig. Qj)J The main effect of increasing M' is the reduc- 
tion of the edge state coupling through the NI strip with 
a corresponding decrease of the k x = minigap at E = 0. 
If we instead decrease M' the minigap at k x — increases 
and the ID edge states are restricted to a smaller spectral 
region since \E\ < M'. 

As a final point, we note that the system with three 
interfaces, having an odd number of helical edge modes 
per spin is protected by TRS from opening a full gap in 
the edge state spectrum notwithstanding the finite over- 



lap of individual edge states. At any given energy within 
the bulk-gap and for any given spin, there is an odd num- 
ber (one or three) of propagating edge modes per spin. 
This is a direct consequence of the conservation of the 
parity of helical edge states linked to the Z2 topolog- 
ical invariant. To put these statements in relation to 
the configurations treated in this work, we note that the 
configuration shown in Fig. QJa) with edge dispersions in 
Fig. 2(a) and the configuration shown in Fig. [JJd) with 
edge dispersions shown in Figs. [12] [HI EU have both an 
odd number of Dirac mass domains (one and three, re- 
spectively) and correspondingly they are metallic. The 
configuration in Fig. HJb) with dispersions in Figs. [5] and 
[9] and the configuration in Fig. QJc) with dispersions in 
Fig. Illf a) have an even number of Dirac mass domains 
and are gapped — correspondingly they are insulators. 



VII. CONCLUSION 

We have analyzed the edge states of a system described 
by the BHZ model where the Dirac mass varies spatially 
thus forming Dirac mass domains where topological in- 
sulating regions alternate with normal insulating regions. 
While for a Tl/vacuum interface the use of hard- wall 
boundary conditions can be assumed, we show that at 
a TI/NI mass domain with a finite NI mass, soft-wall 
boundary conditions (characterized by the continuity of 
the spinor and its derivative) are required to correctly 
account for the edge state dispersion curves and for the 
shape of the corresponding wave functions. We solve the 
edge states for a system up to three TI/NI interfaces. 
For the case of two interfaces, we solve the problem of 
a TI strip with hard-wall, soft-wall and hybrid bound- 
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ary conditions, extending the work in Ref. |3l|. We also 
have investigated the case of edge states that are lat- 
erally coupled via a narrow NI domain. While the edge 
mode spectrum is fully gapped in the two-boundary cases 
due to the edge mode overlap, we show that, as required 
by time-reversal symmetry, in the three boundary sys- 
tem, an odd number of edge modes (one or three) per 
spin is always present at any given energy within the 
bulk gap. The models solved in this work should be rele- 
vant to understand multiple Dirac mass domains induced 
by fractions of nanometer ranged thickness fluctuations 
in HgTe-based quantum wells or via tunable voltage- 
induced band-inversions, e.g. in double quantum well 
structures^ 5 -. Such a controlled creation of multiple heli- 
cal edge states within a single structure could be used to 
create tunable spin- and charge-transport devices^—. 
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Appendix A: The 4-band Model 



If Eq. KA2J) is satisfied but Eq. flMj) is not for 
^min(^min) tnen s * m tne conduction (valence) band has 
a single minimum at k = with energy \M\ (— \M\). 

With standard TI QW parameters the condition in 
Eq. (|A2[) is generally not fulfilled; with our choice of pa- 
rameters it would correspond to 



A 2 

M < — fa -32 meV, 
AB 



(A6) 



and therefore the valence (conduction) band has a max- 
imum (minimum) at k — and energy E = M (E = 
-M). 

Compatible with a fixed energy E and a real k x , one 
generally obtains four complex values of k y , given by 
Eq. ([3]) . Let us now analyze the domain of k y as a func- 
tion of the TI parameters. Eq. §5$ leads to complex so- 
lutions if 



F 2 - Q 2 < 0, 



(A7) 



otherwise the solutions are either purely real or purely 
imaginary. Such an analysis shows that complex fcy's are 
found for E e (e < l,e+) if \B\ > \D\ and in E $ (e^e?) 
for LB I < ID I . 



We consider the usual TI Hamiltonian for one of the 
Kramers partners, given by the following 

h(k)= (C-Dk 2 )a +(M-Bk 2 )<7 z +A{k x a x -k y c7 y ). (Al) 

Let us define the following two conditions 



Note that for \B\ > \D\, if both conduction and valence 
band have a Mexican hat form, then complex fc y 's are 
essentially found inside the gap region, bound by e± in 



Eq. (|A5|) . In this case, ky — k± [defined in Eq. ([3])] are 
both complex in the gap spectral range with \k + \ = \k—\. 
We define 



A 2 
4B' 

{kt m f > 0, 



(A2) 
(A3) 



where v — sg n \B /(D 2 — B 2 )] . If both conditions in 
Eqs. (IA2|) and (|A3[) are satisfied, the bulk energy disper- 
sion, given by Eq. ([2]), has a "Mexican hat" form with a 
local maximum at k — 0, and 



, 2 = M A 2 \A\D / A 2 - 4MB 

B 2B 2 2B 2 V D 2 ~ B 2 ' 1 ' 



are the valence band maxima (— ) and conduction band 
minima (+), respectively, with energies 



4- 



DM DA 2 \A\ — — — , 

X+^±^#MW^). (A5) 



k± = ±k± = ±\]-k 2 x — F± iy/\F 2 -Q 2 \= ±u + iv,(A8) 



where u and v are the real and imaginary parts of k± and 
v > 0. We choose the present definition of k±, so that 
they exponentially decay along the y-axis. Components 
corresponding to k y = k± contribute to edge states (if 
existent) with a single decay length 1/v and an oscillatory 
behavior as sin (uy) . If Eq. (|A2[) is satisfied but Eq. (|A3[) 
is not for t„(fc~ in ) then in the interval between (e^) 
and —\M\ (\M\) k y is purely imaginary. In this case we, 
instead, define k± = k± and, when existent, edge states 
will have two decaying lengths (l/k±) and no oscillatory 
behavio r 15 ' 29 . 
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